Translation of the first example in the talk "Emily Gorcenski - Polynomial Chaos: A technique for modeling uncertainty". The idea is to approximate a general, potentially empirical, probability distribution using a suitable polynomial expansion (Polynomial/Wiener Chaos). With Julia and PolyChaos package it simplified beautifully. For full explanation see related blog post.
#]add PolyChaos, Distributions, Plots
using PolyChaos, Distributions, Plots
First define our Hermite polynomials, $H_i(x)$ (called Gauss in PolyChaos because Hermite is reserved for the physicists Hermite polynomials which are slightly different).
op_gauss = GaussOrthoPoly(20)
H = [x -> evaluate(i, x, op_gauss) for i in 0:5]
plot(-1.5:0.01:1.5, H, leg=false)
inv_cdf(dist) = u -> quantile(dist, u)
h = inv_cdf(Exponential())
l = inv_cdf(Normal())
#7 (generic function with 1 method)
PolyChaos calculates our inner (scalar) product for us, conveniently.
The integrand is just the Galerkin projection with the transformed functions, $h(u)$ and $l(u)$.
sp = computeSP2(op_gauss)
integrand(i) = u -> h(u)*H[i](l(u))
integrand (generic function with 1 method)
Perform the integration up to $p$ polynomials, we can use PolyChaos for this too although any integration scheme works.
p = 5
int_op = Uniform01OrthoPoly(1000, addQuadrature=true)
ki = [integrate(integrand(i), int_op) / sp[i] for i in 1:p]
5-element Array{Float64,1}: 0.9999993692484951 0.9031921803407967 0.29780023616476886 0.03334230632910607 -0.0024169819236717414
Now we have the $k_i$s to do the transform we can reconstitute the approximated distribution by adding up a load of transformed Gaussian random variables, $\zeta_i$.
ζ = randn(5000)
Σ = sum(ki[i] * H[i](ζ) for i in 1:p)
histogram(Σ, normed=true, leg=false)